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Programmed cell death (PCD) is an essential cellular mechanism that is evolutionary con- 
served, mediated through various pathways and acts by integrating different stimuli. Many 
diseases such as neurodegenerative diseases and cancers are found to be caused by, or 
associated with, regulations in the cell death pathways. Yeast Saccharomyces cerevisiae, is 
a unicellular eukaryotic organism that shares with human cells components and pathways 
of the PCD and is therefore used as a model organism. Boolean modeling is becoming 
promising approach to capture qualitative behavior and describe essential properties of 
such complex networks. Here we present large literature-based and to our knowledge first 
Boolean model that combines pathways leading to apoptosis (a type of PCD) in yeast. Analy- 
sis of the yeast model confirmed experimental findings of anti-apoptotic role of Bir1 p and 
pro-apoptotic role of Stm1 p and revealed activation of the stress protein kinase Hog propos- 
ing the maximal level of activation upon heat stress. In addition we extended the yeast 
model and created an in silico humanized yeast in which human pro- and anti-apoptotic reg- 
ulators Bcl-2 family and Valosin-contain protein (VCP) are included in the model. We showed 
that accumulation of Bax in silico humanized yeast shows apoptotic markers and that VCP 
is essential target of Akt Signaling. The presented Boolean model provides comprehen- 
sive description of yeast apoptosis network behavior. Extended model of humanized yeast 
gives new insights of how complex human disease like neurodegeneration can initially be 
tested. 
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INTRODUCTION 

Apoptosis is a complex process which is strictly under control of 
several regulatory networks. Any kind of malfunctioning in these 
controlling systems due to insufficient or excessive apoptosis sig- 
nal can potentially lead to threatening diseases such as various 
types of cancer and neurodegenerative disorders. Therefore keep- 
ing this process tightly regulated is important for the cell. Even 
though apoptosis is often studied in multicellular organisms, the 
discovery of yeast apoptosis in 1997 (Madeo et al., 1997) attracted 
the attention of the wide research community (Frohlich et al., 
2007; Owsianowski et al., 2008; Madeo et al, 2009; Carmona- 
Gutierrez et al., 2010b). As in other multicellular organisms the 
apoptosis in yeast is triggered by both internal and external sig- 
nals. In yeast, the external signals can include acetic acid (Ludovico 
et al., 2001, 2002), salts, metal ions, ethanol, osmotic stress, heat 
stress (Madeo et al, 2009), lipids (Aerts et al, 2008; Low et al, 
2008; Garbarino et al., 2009), mating pheromone (Zhang et al., 
2006), different pharmacological molecules, and drugs (Almeida 
et al, 2008). Internal signals can include, ammonium, NO, ROS 
(that can be generated within the cell by mitochondria and the 
ER, and also induced by H 2 0 2 addition (Madeo et al, 1999) 
and other factors), damage (proteins, lipids, nucleic acids) as a 
consequence of aging and mutations (Mazzoni et al., 2005; Wein- 
berger et al, 2005; Hauptmann et al., 2006), as well as expression 



of heterologous proteins, such as human pro-apoptotic proteins 
(Eisenberg et al., 2007). Many proteins residing in the cytoplasm, 
nucleus, mitochondria, ER, peroxisomes, and lysosomes have been 
identified as the regulators of apoptosis. For example, proteoly- 
sis is one of the main steps that leads to execution of cell death 
and, a yeast metacaspase Ycalp has been shown to be central 
for most (but not all) cell death scenarios (Madeo et al., 2002, 
2009). Besides degradation of proteins, degradation of nucleic 
acids is also carried out during apoptotic death and one of the two 
important caspase-independent mediators is Nuclp (homolog 
of endonuclease G; Buttner et al., 2007). The second is Aiflp 
(apoptosis-inducing factor; Wissing et al., 2004) that is together 
with Nuclp released from the mitochondrion and translocated 
to the nucleus during the initiation and execution of apopto- 
sis. The regulation of apoptosis in the nucleus, is achieved via 
pro-apoptotic factor Nmalllp (nuclear mediator of apoptosis; 
Fahrenkrog et al., 2004) a serine protease that cleaves an anti- 
apoptotic factor (inhibitor of apoptosis, IAP) Birlp, which is 
the only known IAP in yeast, and its anti-apoptotic mechanisms 
(known to be YCA1 -independent) are not well characterized (Wal- 
ter et al, 2006). To understand how large and complex network 
of apoptosis process is regulated it should be studied as a whole 
allowing identification of the properties essential for biological 
function (Janes et al., 2005). To complement experimental studies 
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mathematical models are often use permitting systematic analysis 
of the network components either individually or jointly ( Wolken- 
hauer, 2002; Stelling, 2004). In Boolean networks (BN) introduced 
by Kauffman (1969) these assumptions are made based on activa- 
tion/inhibition effects of one node on another downstream node. 
The Boolean "on" state (or 1 state or "true" state) can be translate 
to biological active state of specific species, while "off" state (or 
0 state or "false" state) corresponds to inactive state. With sim- 
ple logical rules (AND, OR, and NOT) it is possible to capture 
system's behavior in a discrete manner without being dependent 
on experimental measurements such as molecular concentration 
or kinetic rates. This type of model implementation is becom- 
ing more common in biology (Handorf and Klipp, 2012) and 
examples include various model organisms and processes, rang- 
ing from cell cycle models of simple fission yeast (Davidich and 
Bornholdt, 2008), to complex dynamic analysis of mammalian 
cell cycle (Faure et al, 2006), study of mammalian neurotrans- 
mitter signaling pathway (Gupta et al., 2007), investigation of 
irreversible mammalian apoptosis and stable surviving (Mai and 
Liu, 2009), and model of apoptosis in human (Schlatter et al., 
2009). In each of these studies different extensions of Boolean 
modeling is implemented giving clear indication of the grow- 
ing application of logic based modeling in qualitative studies of 
biological networks where there is not much quantitative data 
available. 

We describe here a model based on Boolean network approach 
that consists of two parts: in the first part construction and 
the evaluation of the yeast apoptosis Boolean model is intro- 
duced and in the second part we propose the use of established 
model for study of human apoptotic proteins in yeast. Simplic- 
ity of BN allowed us to construct the model that integrates vast 
amount of heterogeneous knowledge that currently exists for yeast 
apoptosis. The major purpose of this study is to understand the 
emergence of systems properties. Extensive analysis of the state 
space in combination of different input signals generated series 
of in silico experiments. Simulating knock out experiments we 
were able to test the function of specific feedback structures in 
apoptotic network. The results were compared with the exist- 
ing experimental data and the model was used to explore several 
hypotheses in order to better understand certain apoptotic mech- 
anisms and to suggest new strategies for further experimental 
studies. 

RESULTS 

NETWORK TOPOLOGY 

The constructed yeast apoptosis network contains 73 species 
and 115 reactions (Table 1). Six species do not have succes- 
sor (Sink/Output) or predecessor (Source/Input; Table 2). In the 
schematic diagram (Figure 1) species are represented as nodes and 
reactions as edges. Species include: processes, proteins, and small 
molecules (metabolites or signals) and Reactions include activa- 
tion (green arrow) and inhibition (red arrow). Nine nodes depict 
inputs to the model and are colored blue. Seventeen elements are 
active in nucleus and are shown in gray, 12 mitochondria species 
are shown in yellow, and 34 orange boxes represent species resid- 
ing in cytoplasm. System has only one output node which is called 



"Apoptosis" and it is in dark blue. Filled blue circles are used as 
an "AND" gates between two or more reactions to indicate the 
necessity of presence of two species for activation or inhibition of 
a reaction (Figure 1). 

Network is activated via input signals corresponding to the fol- 
lowing species: Acetic Acid, Heat, H2O2, Adozelesin, Mg 2+ , Cu 2+ , 
Salt (NaCl), mating, and osmotic stress. 

MODEL PROPERTIES 

Like many molecular mechanisms in living cells, apoptosis can be 
approximated as an outcome of sequential regulation steps which 
do not occur all at the same time. As a result, upon induction, 
the state of the cell changes with passing time. In order to cap- 
ture this feature, each reaction was assigned activation time scale 
implicating occurrences of different scenarios in sequential order. 
By introducing time delays to the logical function, it is possible 
to describe dynamic behavior of the given process using logical 
networks (Thomas and D'Ari, 1990). 

For this study we used Cell Net Analyser (CNA; Klamt et al, 
2006) which provides the function to capture signal propagation 
in a time series to get a snapshot of the network and discrim- 
ination of signaling events. This approach has been successfully 
used to model human apoptotic network revealing new feedback 
loops (Schlatter et al., 2009). To describe dynamic behavior of 
the yeast apoptosis five timescales t = [0, 2, 4, 5, 6] are assigned. 
These timescales are constants and indicate at which timescale 
each node gets activated or inhibited. Simulation of the network 
at t = x comprises of all interactions with timescale t < x but not 
interactions with t > x. This gives the possibility of separation of 
different functional groups such as different signaling pathways 
or feedback loops. It should be noted that specified timescales do 
not refer to real time but are only indicators of sequential regu- 
latory steps (difference between timescale f = 0 and t = 2 is equal 
to difference between timescales t = 2 and t = 4 or t = 4 and t = 5 
or between t = 5 and t = 6). Timescales are not assigned based 
on speed of reactions (how fast or slow reactions are), they are 
assigned based on the sequence of events. To allow flexibility in 
our model and facilitate changes and insertions of new events 
timescales 1 and 3 are reserved and not used in the current net- 
work structure. First timescale f = 0 is assigned to genes which 
are already present and constantly active in cell (27 reactions). 
Values for all input arcs to stimuli species are set to 1 at second 
timescale t = 2 (53 reactions). Further interactions are activated 
at time t = 4 (17 reactions). Although yeast apoptosis does not 
contain any feedback loops time point t = 5 is reserved for feed- 
back loops which get activated in response to osmotic and heats 
shock and will lead to apoptosis. Considering impact of feed- 
back loops on system behavior it is reasonable to assign them a 
separate timescale. Finally interaction occurring at the very end 
are assigned timescale t = 6 (18 reactions). Details on each reac- 
tion, their time points and species involved in reactions are in 
Table 3. 

The model simulated the induction of apoptosis both in the 
independent mode (assessing each individual stimulus separately) 
as well as in the additive mode, where the activation of all inputs 
were set at the same time (Table 4). 
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Table 1 | List of species in yeast Boolean model. 



ID 


Species 


Type 


Name description 


1 


ABNORMALTELOMERASE 


Change 


ABNORMALTELOMERASE 


2 


ACETIC ACID 


Input 


Chemical 


3 


ADENYLATECYCLASE 


Enzyme 


Lyase enzyme 


4 


ADOZELESIN 


Input 


Drug 


5 


AIF-MT 


Protein 


Apoptosis-inducing factor in mitochondria 


6 


AIF-NUC 


Protein 


Apoptosis-inducing factor in nuc eus 


7 


APOPTOSIS 


Output 


Cell death 


8 


BIR1 


Protein 


Bacu ovira IAP repeat-containing protein 1 


9 


CAMP 


Protein 


Cyclic adenosine monophosphate 


10 


CDC48 


Protein 


Cell division cycle 


11 


CDC6 


Protein 


Cell division cycle 


12 


CPR3 


Protein 


Cyc osporin-sensitive proline rotamase 


13 


CU2 


Input 


Ion 


14 


DesCYCLINCCDK8 


Change 


Destruction of cylinC/CDCk8 


15 


CYTC-CYT 


Protein 


Cytochrome in cytososl 


16 


CYTC-MT 


Protein 


Cytochrome in mitochondria 


17 


DNA-FRAG 


Change 


DNA fragmentation 


18 


DRE2/TAH18 


Change 


Dre2-TAH18 complex 


19 


EMC4 


Protein 


ER membrane protein complex 


20 


ESP1 


Protein 


Separase 


21 


FIS1 


Protein 


Mitochondrial FISsion 


22 


FYV10 


Protein 


Function required for yeast viability 


23 


H2B 


Protein 


Histon 2B 


24 


H202 


Input 


Hydroxide peroxide 


25 


HEAT 


Input 


Event 


26 


HK 




House keeping function** 


27 


HOG1 


Protein 


High osmo arity glycerol response 


28 


HOG1-DEP 


Protein 


HOG1 dependent genes 


29 


HOS3 


Protein 


Hda one similar 


30 


KAP123 


Protein 


KAryoPherin 


31 


MAPK 


Protein 


Map kinase cascade 


32 


MATING 


Input 


Mating pheromone 


33 


MCD1-MT 


Protein 


Mitotic chromosome determinant in mitochondria 


34 


MCD1-NUC 


Protein 


Mitotic chromosome determinant in nuc eus 


35 


MDV1 


Protein 


Mitochondrial Division 


36 


MEC1 


Protein 


Mitosis entry checkpoint 


37 


MG2 


Input 


Ion 


38 


MMI1 


Protein 


Translation machinery associated 


39 


MSN2-4 


Protein 


Mu ticopy suppressor of SNF1 mutation 


40 


MT-ALT 


Change 


Mitochondria alteration 


41 


MT-FRAG 


Change 


Mitochondria fragmentation 


42 


NDI1 


Protein 


NADH dehydrogenase internal 


43 


NMA111-CYT 


Protein 


Nuclear mediator of apoptosis in cytososl 


44 


NMA111-NUC 


Protein 


Nuclear mediator of apoptosis in nucleus 


45 


NUC1-MT 


Protein 


NUCIease 1 in mitochondria 


46 


NUC1-NUC 


Protein 


NUCIease 1 in nucleus 


47 


PKA 


Protein 


Protein kinase A 


48 


POR1-2 


Protein 


PORin 


49 


PROTOSOM 


Complex 


PROTOSOM 


50 


PTP2 


Protein 


Protein tyrosine phosphatase 



(Continued) 
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Table 1 | Continued 



ID 


Species 


Type 


Name description 


51 


PTP3 


Protein 


Protein tyrosine phosphatase 


52 


RAS2 


Protein 


Homologous to RAS proto-oncogene 


53 


RedActinDyn 


Change 


Reduced actin dynamic 


54 


RLM1 


Protein 


Resistance to lethality of MKK1 P386 overexpression 


55 


ROS-OY I 


Molecule 


Reactive oxygen species in cytososl 


56 


ROS-MT 


Molecule 


Reactive oxygen species in mitochondria 


57 


RPD3 


Protein 


Reduced potassium dependency 


58 


SALT 


Input 




59 


SDP1 


Protein 


Stress-inducible dual specificity phosphatase 


60 


SLT2 


Protein 


SYNtaxin (SYN8) 


61 


SN01 


Protein 


SNZ proximal open reading frame 


62 


SOD1 


Protein 


Superoxide dismutase 


63 


SOD2 


Protein 


Superoxide dismutase 


64 


SR07 


Protein 


Suppressor of rho3 


65 


STE20-CYT 


Protein 


Sterile in cytosol 


66 


STE20-NUC 


Protein 


Sterile in nucleus 


67 


STM1-CYT 


Protein 


Translation initiation factor (TIF3) in cytosos 


68 


STM1-NUC 


Protein 


Translation initiation factor (TIF3) in nucleus 


69 


STRESS 


Input 


Event 


70 


SVF1 


Protein 


Survival factor 


71 


TAT-D 


Protein 


3'^ 5' exonuclease and endonuclease 


72 


TOR1 


Protein 


Target of rapamycin 


73 


YCA1 


Protein 


MetaCAspase 



Table includes name, type, and description of each species involved in yeast apoptosis. 

'Species refers to all types of nodes that are depicted in the network map and are color coded based on their location or activity in cell. 
'"Housekeeping function refers to those genes which are present and give snapshot of state of the cell before applying any kind of treatment. 



of apoptosis by activation of Hoglp by heat stress, (ii) inhibition 
by Birlp of acetic acid-induced apoptosis, and (iii) induction of 
apoptosis with H2O2 by activation of Stmlp. We show here that 
the behavior that emerges from specific interactions in the model 
is in agreement with published experimental data. 

ACTIVATION OF HEAT STRESS-INDUCED APOPTOSIS WITH Hoglp 

We simulated the activation of the mitogen-activated protein 
kinase that is a key component of the HOG pathway (Albertyn 
and Hohmann, 1994; Van Wuytswinkel et al, 2000; Hohmann, 
2002; de Nadal et al., 2004), which is an osmoregulatory sig- 
nal transduction cascade (Hohmann, 2009). Upon stress, Hoglp 
(encoded by HOG1/YLR1 13W) regulates the expression of almost 
600 genes by phosphorylating several different transcription fac- 
tors (Hohmann, 2002; Westfall et al, 2004). The activity and 
nuclear localization of Hoglp is regulated by its phosphorylation 
state, and that in turn is regulated by the kinase MAPKK Pbs2p and 
the phosphatases Ptclp, Ptc2p, Ptc3p, Ptp2p, and Ptp3p (Brewster 
et al, 1993; Ferrigno et al, 1998; Warmka et al., 2001; Young et al, 
2002). Besides the induction by osmotic stress, the HOG path- 
way can be induced by heat stress, via a Sholp-dependent sensory 
mechanism (Winkler et al, 2002), thus we used heat stress as an 
input to activate the HOG pathway and simulate the cell death 
response. Heat stress activates Slt2 which then activates Rlml. 
Consequently Rlml triggers expression of Slt2 and PTP2 form- 
ing a feedback loop. On the other hand heat shock inhibits activity 



Table 2 | Summary of the interactions without successor 
(Sink/Output) and predecessor (Source/Input). 



Species 


Type of connection 


Number of connections 


Ros-MT 


Sink/output 


1 


MCD1-NUC 


Sink/output 


1 


H2B 


Sink/output 


5 


CAMP 


Sink/output 


2 


RedActinDyn 


Source/input 


3 


AbnormalTelomer 


Source/input 


1 


HK 


Source/input 


18 



Table contains species without successor (there are no edges coming out of these 
nodes! and predecessor (there are no edges going to these nodes). 



PREDICTIONS WITH THE CONTINUOUS MODEL 

Based on qualitative knowledge we have constructed a discrete 
Boolean model of yeast apoptosis. While it can capture its essen- 
tial behavior, the question remained how this model can be used to 
predict the qualitative behavior of the system. In order to address 
this question we expanded the model by transforming the discrete 
Boolean model into a continuous model using Odefy (Wittmann 
et al, 2009; Krumsiek et al, 2010; see Materials and Methods). 
In order to test the continuous model and its predictive capac- 
ity, we performed three independent case-studies: (i) induction 
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Adozelesin H202 
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Mating — Heat 





■ Inputs 
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Mitochondria 


■ 


Nucleus 


* Inhibition 



FIGURE 1 | Schematic representation of Yeast Apoptosis network. Blue boxes depict input nodes, yellow nodes are placed in mitochondria, orange nodes 
reside in cytosol and grey nodes belong to nucleus. Green arrows show activation effect and red arrows show inhibition effect. Blue circles depict "AND" gate. 



of PICA leading to release the inhibition effect of PKA on MSN2- 
4. MSN2-4 activates Sdpl which then along with PTP2 inhibit 
Slt2. We started by setting all other stimuli to zero and then by 
transforming the Boolean apoptosis model into Hill Cube contin- 
uous model we performed the simulations. Steady states predicted 
from continuous Hill Cube (Figure 2) and synchronous Boolean 
(Figure 3) model are in perfect agreement with each other for all 
nodes apart from the nodes corresponding to the feedback loops 
in response to heat and osmotic stress. Feedback loop in heat acti- 
vated pathway includes activation of Rlml by Slt2 and expression 
of Slt2 by Rlml. Osmotic shock induced Hogl activates Rlml 
to regulate Slt2 which was inhibited by osmotic shock and PTP2. 
Unknown Hogl dependent transcription factor triggers transcrip- 
tion of PTP3 and PTP2 resulting in dephosphorylation of Hogl 
phosphotyrosine which inhibits Hogl activities. All nodes in our 
model were connected to activation of apoptosis via the heat shock 
pathway and were indeed activated in our simulation and are col- 
ored in red; other pathways are colored in blue indicating that they 
are not activated during the simulation (Figure 2). Model predicts 
that upon heat induction, concentration of Hoglp changes trough 
time never reaching its maximum level of concentration and as an 
intensity of stimulus decreases level of Hoglp also decreases. Our 



simulations also proposed that the maximal level of activation of 
Hoglp during heat stress is 70% of the total activation and as the 
heat stimulus continues over time its activity decreases and reaches 
a plateau at 40% of total activation (Figure 4). 

This model prediction is in agreement with an experimental 
measurement performed by Winkler et al. (2002) which shows the 
induced activity of Hoglp upon heat induction (Figure 4). The 
difference between the maximal activity values and the lag phase of 
the activation is due to the fact that we did not use any quantitative 
data as input to the model. 

INHIBITION OF ACETIC ACID-INDUCED APOPTOSIS WITH Birlp 

Even though there are several genes and pathways that are involved 
in yeast pro-apoptotic response, there are only few anti-apoptotic 
genes. For many years Birlp (encoded by BIR1/YJR089W) was 
thought to be the only apoptotic inhibitor in yeast, but recently 
a TSC22 -protein family was found and two of the proteins 
with TSC22-motif have been shown to have anti-apoptotic roles 
Snolp and FvylOp (Khoury et al, 2008). Birl belongs to the 
IAP family and phylogenetic analysis (Uren et al, 1998) showed 
similarity to Schizosaccharomyces pombe BIR1, human survivin, 
and Caenorhabditis elegance BIR-1 and BIR-2 proteins. Yeast's 
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Table 3 | List of logical interactions in yeast Boolean model. 



ID 


Interaction 


Function 


Time scale 


Reference 


1 


HK = AIF1-MT 


Housekeeping 


0 




2 


HK = DRE2/TAH18 


Housekeeping 


0 




3 


HK= EMC4 


Housekeeping 


0 




4 


HK= SVF1 


Housekeeping 


0 




5 


HK= FVY10 


Housekeeping 


0 




6 


HK = SOD2 


Housekeeping 


0 




7 


HK = SN01 


Housekeeping 


0 




8 


HK= NDI1 


Housekeeping 


0 




9 


HK= POR1-2 


Housekeeping 


0 




10 


HK= MMI1 


Housekeeping 


0 




11 


HK= MCD1-MT 


Housekeeping 


0 




12 


HK=SR07 


Housekeeping 


0 




13 


HK=CDC48 


Housekeeping 


0 




14 


HK= FIS1 


Housekeeping 


0 




15 


HK= MDV1 


Housekeeping 


0 




16 


HK = STM1-CYT 


Housekeeping 


0 




17 


=AceticAcid 


Input 


2 




18 


=Adozelesin 


Input 


2 




19 


=CU2 


Input 


2 




20 


=H 2 0 2 


Input 


2 




21 


=Mating 


Input 


2 




22 


=MG2 


Input 


2 




23 


=Salt 


Input 


2 




24 


=Heat 


Input 


2 




25 


=Stress 


Input 


2 




26 


ISOD2 + NDI1 =ROS-MT 




4 


Li etal. (2006) 


27 


AceticAcid = CytC-MT 




4 


Ludovico et al. (2002) 


28 


CDC48 = CytC-CYT 




4 


Eisenberg et al. (2007) 


29 


CytC-CYT=YCA1 




4 


Eisenberg et al. (2007) 


30 


CytC-MT = CytC-CYT 




4 


Eisenberg et al. (2007) 


31 


MCD1-MT = CytC-MT 




4 


Yang et al. (2008) 


32 


MEC1 =YCA1 




4 


Weinberger et al. (2005) 


33 


MT-Frag = MT-ALT 




4 


Wissing etal. (2004) 


34 


I FYV10 = Apoptosis 




4 


Khoury etal. (2008) 


35 


2 CDC48 = ROS-CYT 




4 




36 


CU2 + CPR3= Apoptosis 




4 


: : and Zhou (2007) 


37 


DNA-Frag = Apoptosis 




4 


Madeo et al. (2009) 


38 


ESP1 = ROS-CYT 




4 


Yang et al. (2008) 


39 


MT-Frag =YCA1 




4 


Eisenberg et al. (2007) 


40 


NMA111-CYT= NMA111-NUC 




4 


Walter etal. (2006) 


41 


NUC1-MT = KAP123 




4 


Buttneretal. (2007) 


42 


RAS2 = AdenylateCyclase 




4 


Wood et al. (1994) 


43 


RAS2 = ROS-CYT 




4 


Kataoka et al. (1984) 


44 


RedActinDyn = ROS-CYT 




4 


Eisenberg etal. (2007) 


45 


RedActinDyn=YCA1 




4 


Madeo et al. (2009) 


46 


ROS-CYT = Apoptosis 




4 


Eisenberg et al. (2007) 


47 


ROS-CYT =YCA1 




4 


Madeo et al. (2002) 


48 


Salt = ROS-CYT 




4 


Wadskog et al. (2004) 


49 


SOD1 = ROS-CYT 




4 


Eisenberg et al. (2007) 


50 


STE20-NUC=H2B 




4 


Madeo et al. (2009) 



(Continued) 
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Table 3 | Continued 

ID Interaction Function Time scale Reference 



51 Stress = RPD3 4 Ahn et al. (2006) 

52 Tat-D = DNA-Frag 4 Qiu et al. (2005) 

53 ISN01 =Apoptosis 4 Khoury et al. (2008) 

54 AbnormalTelomer = MEC1 4 Weinberger et al. (2005) 

55 AdenylateCyclase = CAMP 4 Schmelzle et al. (2004) 

56 Adozelesin = CDC6 4 Blanchard et al. (2002) 

57 AIF1-MT = AIF1-NUC 4 Wissing et al. (2004) 

58 AIF1-NUC = H2B 4 Wissing et al. (2004) 

59 Apoptosis = Output 4 

60 CDC6 = Protosom 4 Blanchard et al. (2002) 

61 ESP1 = MCD1-NUC 4 Yang et al. (2008) 

62 H 2 0 2 = NUC1-MT 4 Buttner et al. (2007) 

63 H 2 0 2 = ESP1 4 Yang et al. (2008) 

64 H 2 0 2 = HOS3 4 Carmona-Gutierrez et al. (2010a) 

65 Heat= NMA111-CYT 4 Walter et al. (2006) 

66 HOS3 = H2B 4 Carmona-Gutierrez et al. (2010a) 

67 KAP123 = NUC1-NUC 4 Buttner et al. (2007) 

68 MAPK = STE20-CYT 4 Carmona-Gutierrez et al. (2010a) 

69 MatingPheromone = MAPK 4 Carmona-Gutierrez et al. (2010a) 

70 Mg 2 +=Tat-D 4 Qiu et al. (2005) 

71 MMI1 = MT-ALT 4 Eisenberg et al. (2007) 

72 MT-ALT =MT-FRAG 4 

73 NUC1-NUC = H2B 4 Buttner et al. (2007) 

74 PKA = MT-ALT 4 Carmona-Gutierrez et al. (2010a) 

75 RAS2 = MT-ALT 4 Eisenberg et al. (2007) 

76 RAS2 = PKA 4 Carmona-Gutierrez et al. (2010a) 

77 RedActinDyn = RAS2 4 Eisenberg et al. (2007) 

78 RPD3 = H2B 4 Ahn et al. (2006) 

79 STE20-CYT = STE20-NUC 4 Carmona-Gutierrez et al. (2010a) 

80 Stress =AdenylateCyclase 4 Schmelzle et al. (2004) 

81 TOR1=CAMP 4 Schmelzle et al. (2004) 

82 TOR1=RAS2 4 Schmelzle et al. (2004) 

83 Heat = SOD1 4 

84 2 NDI1 =ROS-CYT 4 

85 Stress =TOR1 4 Schmelzle et al. (2004) 

86 !PTP2 = SLT2 5 Hahn and Thiele (2002) 

87 !PTP2 = HOG1 5 Hahn and Thiele (2002) 

88 HOG1 =HOG1-Dep 5 Hahn and Thiele (2002) 

89 !PTP3 = HOG1 5 Hahn and Thiele (2002) 

90 !SDP1=SLT2 5 Hahn and Thiele (2002) 

91 !SLT2 = DesCyclinCCDK8 5 Krasley et al. (2006) 

92 !Stress = SLT2 5 Hahn and Thiele (2002) 

93 DesCyclinCCCDK8 = ROS-CYT 5 Krasley et al. (2006) 

94 Heat=PKA 5 Hahn and Thiele (2002) 

95 Heat = SLT2 5 Hahn and Thiele (2002) 

96 Hog1 = RLM1 5 Hahn and Thiele (2002) 

97 HOG1-Dep = PTP3 5 Hahn and Thiele (2002) 

98 MSN2-4 = SDP1 5 Hahn and Thiele (2002) 

99 PKA = MSN2-4 5 Hahn and Thiele (2002) 

100 RLM1 = PTP2 5 Hahn and Thiele (2002) 

101 RLM1 = SLT2 5 Hahn and Thiele (2002) 



(Continued) 
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Table 3 | Continued 



ID 


Interaction 


Function 


Time scale 


Reference 


102 


SLT2 = RLM1 




5 


Hahn andThiele (2002) 


103 


Stress = HOb1 




5 


Hahn and Ihiele (2002) 


104 


ISR07 + Salt=YCA1 




6 


Wadskog et al. (2004) 


105 


ISTM1-NUC = DNA-Frag 




6 


Ligr et al. (2001) 


106 


AceticAcid + ISVF1 = ROS-CYT 




6 


Vander Heiden et al. (2002) 


107 


H2O2 + I F_MC4 = ROS-CY I 




6 


Ring et al. (2008) 


108 


H2O2 + ISVF1 = ROS-CYT 




6 


Vander Heiden et al. (2002) 


109 


!FIS1 + MDV1 = MT-Frag 




6 


Eisenberg et al. (2007) 


110 


! PO R 1-2 + AceticAcid = Apoptosis 




6 


Pereira etal. (2007) 


111 


!POR1-2 + H 2 0 2 =Apoptosis 




6 


Pereira et al. (2007) 


112 


H 2 0 2 + !DRE2/TAH18 = MT-Frag 




6 


Vernis et al. (2009) 


113 


YCA1 + IBIR1 = Apoptosis 




6 


Walter etal. (2006) 


114 


!NMA111-NUC = BIR1 




6 


Walter et al. (2006) 


115 


STM1-CYT+ IProtosom = STM1-NUC 




6 


Ligr etal. (2001) 



Table includes involved species in each interaction, logical rule for each interaction and time scale in which interaction takes place. Note: the interactions of the model 
are given in the notation of the Cell Net Analyser: a logical NOT is represented by ".' "; a logical AND is represented by "+ " and interaction on right hand side of "= " 
gives the value of node on left hand side or equation. 



Table 4 | Predicted states of relevant species at steady state. 
Species f = 0 t = 4 ( = 5 



f = 6 



Acetic acid 

Adozelesin 

AIF1-MT 

AIF1-NUC 

Apoptosis 

BIR1 

DNA-FRAG 

H202 

Heat 

NMA111-NUC 

ROS-CYT 

STE20-CYT 

STM1-CYT 

STM1-NUC 

YCA1 



Table includes time scale scenario upon induction of all stimuli. T stands for differ- 
ent time steps. As an example if Acetic Acid is inducted to the cell at time point 0 
it is expected to cause DNA Fragmentation at next time step and eventually ends 
to apoptosis in last step. 

Birlp has been previously intensively studied in chromosome seg- 
regation (as a component of the Aurora kinase complex (chromo- 
some passenger complex, CPC; Ruchaud et al., 2007) but recently 
studies of its role as a negative regulator of apoptosis have gained 
momentum. It has been shown that Birlp is a target for degrada- 
tion by Nmalllp (Owsianowski et al., 2008), an apoptotic serine 
protease and yeast cells lacking Birl are more sensitive to apopto- 
sis, while overexpression of Birl reduces apoptosis (Walter et al, 
2006). 



When building the model we assumed that addition of 
acetic acid induces cytochrome c release from mitochondria 
and its translocation to the cytosol. Another assumption is that 
the execution of apoptosis is Ycalp-dependent (encoded by 
YCA1/MCA1/YOR197W) in this context and that the activation 
of the caspase is downstream from the cytochrome c transloca- 
tion (Ludovico et al, 2001, 2002). Thus, taking these assumptions 
into account acetic acid was used as an inducer in the simula- 
tions with direct induction of apoptosis. The model confirmed 
experimental evidence that Birl is indeed apoptosis inhibitor in 
yeast (Walter et al, 2006) and as acetic acid is applied as a pulse 
stimulus (that is then decreasing over time) and cytochrome c 
present in the cytosol increases, the decrease of Birlp due to 
Nmalllp degradation promotes apoptosis that then reaches the 
maximum (Figure 5A). We then performed an artificial and bio- 
logically irrelevant simulation: after apoptosis has occurred and 
no more inducer was present (acetic acid is not added again and 
the previously added amount has been "used") the eventual the- 
oretical accumulation of Birlp did inhibit apoptosis and revert it 
to zero. Obviously this scenario is biologically implausible, since 
once the cell has undergone apoptosis it cannot be revived, but 
has a purpose to show mathematical correctness of the developed 
model. 

Additionally, we have converted discreet apoptosis model to 
continuous model with constant presence of Birl in order to 
test if this type of model would predict the same outcome as 
the experimental approach in which Birlp is overexpressed and 
provides protection from induction of apoptosis. The continu- 
ous model was created taking into account translocation effect 
of each element from one compartment to another. This effect 
has been applied to translocation of Nucl from mitochondria to 
nucleus, MCD1 from nucleus to mitochondria, Stml from cytosol 
to nucleus Aifl (from mitochondria to nucleus, NMA111 from 
cytosol to nucleus, Ste20 from cytosol to nucleus, and cytochrome 
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FIGURE 2 | Continuous Hill Cube transformation. All species in apoptosis network are mapped on vertical axis. Dark blue color indicates those nodes that are 
not activated (value = 0) while dark red refers to nodes which are completely activated (value = 1). Each color in between indicates the level of activation 
between 0 and 1. 



c from mitochondria to cytosol. Initially, in this model, translo- 
cation of each species is modeled as one node in order to rep- 
resent the re-localization which introduced self-regulatory loops 
to the system which is impossible since to our knowledge yeast 
apoptosis regulatory network does not have such loops. One 
solution implies intuitive approach in defining single variable 
for single species to mimic the effect of transferring from one 
compartment to another. Technically this was solved by intro- 
ducing two variables for a single species representing both com- 
partments they can belong to. This approach, also, confirmed 
experimental finding (Walter et al, 2006) that constant pres- 
ence of Birl, inhibits apoptosis, validating Birl anti-apoptotic role 
(Figure 5B). 

ACTIVATION OF H 2 0 2 -INDUCED APOPTOSIS WITH Stmlp 

It has been shown that degradation of short-lived pro-apoptotic 
proteins via proteasomal - ubiquitine pathway plays important 
role in mammalian apoptosis (Drexler, 1997). One of the yeast 
proteasomal substrates - Stml (YLR150W) was identified in the 
study by Ligr et al. (2001) as an activator of the cell death process. 
Conserved orthologs of Stml are detected in several species (highly 



conserved ortholog in Schizosaccharomyces pombe and a putative 
ortholog in Drosophila melanogaster (Nelson et al, 2000) suggest- 
ing that regulation of apoptosis via this protein can be evolutionary 
conserved process. 

We have in silico tested two scenarios confirming the role of 
Stml in apoptosis (Ligr et al, 2001). Accumulation of Stml will 
induce apoptotic behavior followed by DNA fragmentation which 
is a known marker of cell death, while the Stml knock out will 
promote survival and consequently deterring DNA fragmentation 
(Figures 6A,B). 

IN SILICO HUMANIZED YEAST APOPTOSIS 

S. cerevisiae is a model organism that has conserved genes, pro- 
teins and pathways that are similar to the ones in the human cells. 
This allows for using yeast as a host in which human genes can be 
expressed, so called "humanized yeast" (Winderickx et al., 2008) 
and subsequently the physiological roles and molecular mecha- 
nisms can be studied. In order to test if our system could be used 
in a similar way (in silico humanized yeast) and provide reliable 
simulations and predictions, we inserted three human genes with 
complete downstream pathways into our initial yeast apoptosis 
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FIGURE 3 | Synchronous Boolean model. All species in apoptosis network are mapped on vertical axis. Dark blue indicates those nodes that are not activated 
(value = 0) while dark red refers to nodes which are totally activated (value = 1). 



Boolean model. Yeast apoptotic network was "humanized" by 
in silico insertion of genes belonging to Bcl-2 protein family 
and Valosin-contain protein (VCP) - distinctive representatives 
of human apoptosis (Figure 7). 

INSERTION OF Bcl-2 FAMILY PROTEINS: ANTI-APOPTOTIC Bcl-xL AND 
PRO-APOTOTIC BAX 

This protein family consists of apoptotic agonist and anti-agonist 
proteins. Among members of Bcl-2 family Bax as an apoptotic 
inducer and Bcl-xL as an anti-apoptotic factor were included in 
our extended model. Bax is located in cytosol in mammalian cells 
and has vital role in mitochondria morphogenesis. Heterologously 
expressed Bax causes growth arrest and rapid cell death in S. cere- 
visiae (Sato et al., 1994; Greenhalf et al., 1996). Also expression 
of Bax has been linked to the release of cytochrome c from mito- 
chondria (Manon et al, 1997). Yeast cells with Bax expression 
accumulate ROS and show other apoptotic hallmarks such as DNA 
fragmentation (Madeo et al, 1999). 

Since members of this protein family were previously success- 
fully expressed in yeast, we first validate the "humanized" yeast 
model by inserting Bcl-xL, Bax, Bad, and Bcl-2 genes (Sato et al., 
1994; Greenhalf et al., 1996; Madeo et al., 1999). As an input signal 



human Akt signaling was used (Akt Signaling), as an extrinsic 
regulatory switch (since this signaling cascade is not present in 
yeast). The Akt Serine/Threonine-kinase promotes cell survival by 
phosphorylating the pro-apoptotic protein BAD (member of the 
Bcl-2 family), which is the cause of dissociation of BAD from the 
Bcl-2/Bcl-X complex, and promotion of cell survival (Datta et al, 
1997). Besides cell survival, Akt signaling is related to the cell cycle, 
metabolism, and angiogenesis and therefore a target for anticancer 
drug development (Falasca, 2010; Hers et al, 2011). 

Upon activation of Akt Signaling at the first time step (t = 0, 
node value Akt Signaling = 1), Bcl-2 and Bcl-xL are activated in 
the following step (t = 4, node value Bcl-xL = 1 and Bcl-2 = 1) 
and remain active throughout the simulation and as the inhibitors 
of apoptosis promote the survival (node value Apoptosis = 0; 
Table 5). Simulation results suggest the anti-apoptotic role of two 
Bcl-2 family members: Bcl-2 and Bcl-xL which is consistent with 
experimental findings (Kharbanda et al., 1997). 

DUAL FUNCTIONALITY OF VCP IN SURVIVAL AND APOPTOSIS 

The evolutionary conserved Valosin-containing protein (VCP) is 
a mammalian ortholog of yeast Cdc48, which is the first apoptotic 
mediator found in S. cereivisae (Braun and Zischka, 2008). VCP is 
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the member AAA-ATPase family which is ubiquitously expressed 
(Braun and Zischka, 2008). It consists of four domains: N-terminal 
domain, two ATPase domains Dl and D2 and C-terminal domain 
(Wang et al., 2004). The major ATPase activity of VCP is carried 
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FIGURE 4 | Hog1 study. Comparison of experimental and simulation study. 
Black curve shows activity of wild type Hog1 upon induction of heat 
(Winkler et al., 2002) and red curve illustrates Hog1 activation level in 
continues simulation. 



on by D2 domain (Wang et al, 2004; Song et al, 2003). VCP is 
involved in different cellular process such as protein degradation, 
membrane fusion and chaperone activity (Braun and Zischka, 
2008). Role of VCP/Cdc48 in fluctuating number of death cell 
in various types of disease naming cancer and protein deposition 
diseases is not well understood. Increase in expression of VCP 
is correlated to the development of cancer and metastasis there- 
fore detecting the level of VCP expression is proposed as cancer 
progression marker. Also VCP is known as detector of aggregated 
proteins causing neurodegenerative disease such as Parkinson and 
Alzheimer (Hirabayashi et al., 2001; Mizuno et al., 2003; Ishigaki 
et al., 2004). It has been observed that yeast cell carrying mutation 
of Cdc48 is showing morphological markers of apoptosis (Madeo 
et al., 1997) and that depletion of this protein has apoptotic effect 
in other organisms as well (Imamuraab et al, 2003). 

It had been shown that VCP can both promote or inhibit apop- 
tosis (Braun and Zischka, 2008). Deletion of VCP triggers ER stress 
which is followed by unfolded protein response and cell death via 
ER-associated degradation (ERAD) pathway - a highly conserved 
pathway between mammalian and yeast cells (Braun and Zischka, 
2008). Wild type VCP has pro-apoptotic role in the cells undergo- 
ing apoptosis upon ER stress which in turn triggers ER-associated 
apoptotic pathway. Finally VCP can trigger survival pathway in 
response to NFkB which is a pro-survival molecule in cells with 
over expressed level of VCP/Cdc48(Braun and Zischka, 2008). 

The dual functionality of VCP in survival and apoptosis mech- 
anism makes it an attractive candidate to be inserted in initial 
yeast apoptosis Boolean network. Since VCP includes expression 
of a handful of heterologous proteins, its in vivo insertion to 
yeast plasmid would be difficult to perform and thus represents a 
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FIGURE 5 | Birl study. (A) Acetic acid (green) is applied as a pulse stimulus (that is then decreasing over time) and cytochrome c (yellow) present in the cytosol 
increases, the decrease of Birl p (red) promotes apoptosis (B) Constant presence of Birl p (red) inhibits apoptosis (black), validating Birl anti-apoptotic role. 
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FIGURE 6 | Stml study. (A) Presence of Stmlp (blue) promotes apoptosis (yellow; B) Knock out of Stmlp (blue) prevents apoptosis (yellow) and DNA 
fragmentation (red) and consequently promotes survival. 
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FIGURE 7 | In silico "humanized yeast apoptosis network" Human apoptotic pathways BCL-2 protein family and VCP dependent genes inserted into the 
yeast apoptosis network (due to simplicity we show only human pathways). 
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Table 5 | Summary of simulations upon insertion of Bcl-2 pathway. 



Species 


f = 0 


( = 4 


t = 5 


f = 6 


UV 


0 


0 


0 


0 


Akt signaling 


1 


1 


1 


1 


BCL-XL 


0 


1 


1 


1 


P53 


0 


0 


0 


0 


BAD 


0 


0 


0 


0 


BAX 


0 


0 


0 


0 


BCL-2 


0 


1 


1 


1 


Apoptosis 


0 


0 


0 


0 



Table Includes simulation results of heterologous expression of BCL-2 in yeast 
apoptosis model. T stands for different time steps. Upon activation of Akt Signal- 
ing pathway at first time step BCL-2 gets activated in following step and as an 
inhibitor of apoptosis prevents apoptosis till end point. 



good candidate for further investigation and exploration of model 
capabilities. 

In our model, yeast Cdc48 was replaced with human VCP gene 
and its downstream pathway effector caspases (caspases 3, 6 and 
7), initiator caspases (caspase9 and 12), IAP family, IkB-alpha 
inhibitory protein which when degraded by proteasome cause the 
release, and nuclear translocation of active NFkB (represented in 
the model as NFkB_Cyt and NFkB_Nuc) and gp 130/Stat3 pathway 
(represented as a single node). 

Upon activation of Akt Signaling pro-survival role of VCP is 
observed (node Survival = 1 at t = 6, Table 3). This is achieved 
when IkBa is dissociated from NFkB (node NFkB/IkBa = 1 at 
f = 4), which in turn gets activated (node NFkB=l, f = 5) 
and is translocated to the nucleus (nodes NFkB_Cyt = 1 and 
NFkB_Nuc=l, t = 6) inducing survival of the cell (node Sur- 
vival = 1, t = 6). Simultaneously, activation of NFkB, promotes 
activity of IAP family of proteins which inhibits the activity of both 
effector and initiator caspases (nodes c-3-6-7 = 0 and c-9-12 = 0, 
t = [0,4,5,6] , thus disabling apoptosis-dependent caspase pathway. 
As an independent, but parallel process, activation of VCP via 
Akt Signaling activates gpl30-Stat3 (node gpl30-Stat3 = 1, t = 4) 
pathway and consequently leads to survival (Table 6). 

This result suggests that yeast carrying mammalian VCP gene 
exhibits the same behavior as it is known from the human 
apoptotic model (Vandermoere et al., 2006). 

DISCUSSION 

In this work we have constructed a Boolean model for the bio- 
chemical network that controls apoptosis pathway in budding 
yeast Saccharomyces cerevisiae. Firstly, we presented the Boolean 
model describing only yeast pro and anti-apoptotic genes and val- 
idated the model by further analyzing the role of Stmlp, Birlp, 
and Hoglp. Even though construction of the yeast apoptosis net- 
work involved certain simplifications (nodes have only two states 
and rules are describing network dynamics) we were able to model 
more complex network than using dynamic modeling approach 
that requires knowledge of kinetic parameters for all molecular 
processes. This approach was able to suggest general design prin- 
cipals of yeast apoptosis. One of the advantages of constructing 
mathematical models in biology is that we are able to simulate 



Table 6 | Summary of simulations upon insertion of VCP pathway. 



Species 


t = 0 


t = 4 


f = 5 


f = 6 


Akt signaling 


1 


1 


1 


1 


GP130-STAT3 


o 


1 


1 


1 


NFkB-CYT 


o 


1 


1 


1 


NFkB-NUC 


0 


1 


1 


1 


NFkB/IkBa 


0 


1 


1 


1 


VCP 


1 


1 


1 


1 


C-3-6-7 


0 


0 


0 


0 


C-9-12 


0 


0 


0 


0 


IkBa 


0 


1 


1 


1 


IAP 


0 


1 


1 


1 


Survival 


0 


0 


0 


1 


Apoptosis 


0 


0 


0 


0 



Table includes simulation results of humanized yeast apoptosis model by inser- 
tion of VCR T stands for different time steps. VCP in presence of Akt Signaling is 
expected to prevent apoptosis and promote survival at last time point. 



scenarios that are not feasible in real experiments. In the case 
of Birl p we were able to in silico revive a cell already reaching 
apoptosis. These examples are important to understand general 
mechanisms of the constructed network. A study of proteasomal 
substrate Stmlp suggested that Stmlp seems to be an effector of 
H2O2 and the presence of Stmlp with an inducer leads to apopto- 
sis in our model, and in Stmlp knock out survival is promoted and 
DNA fragmentation is avoided, irrespective of the presence of the 
inducer. Despite the limitation of BN in terms of giving quanti- 
tative predictions of the system dynamics they allow investigation 
of the large networks and their systematic exploration resulting 
in better understanding of cellular processes. Example of Hoglp 
study showed that our model was capable of reproducing the pat- 
tern of Hoglp activation, but was not able to quantitatively predict 
the maximal Hoglp activity. 

In the second stage, we extended the initial model by inserting 
two pathways of human apoptosis (VCP and Bcl-2 family), hereby 
creating a "humanized" in silico yeast. Humanized yeast strains are 
used in experimental molecular biology where the human proteins 
(causing, or associated with diseases) are expressed and studied 
in vivo (in yeast in this case). 

Conservation of many apoptotic mediators and mechanisms 
among the Eukarya provides the possibility of inserting genes 
from other organisms into yeast (in vivo) or in the yeast apop- 
totic network (in silico). We validated constructed humanized 
yeast model by insertion of Bcl-2 family which have been pre- 
viously successfully expressed in yeast. The model was able to 
reproduce well-known pro and anti-apoptotic phenotypes con- 
firming that yeast expressing Bax accumulated ROS and showed 
other apoptotic markers like DNA fragmentation. To test whether 
or not specific pathways playing role in neurodegeneration and 
cancer can be elucidated by our Boolean model and can we con- 
sequently take advantage of a simple system like yeast to explore 
hypothesis generated for higher organisms, we in silico expressed 
evolutionary conserved VCP and its downstream components in 
existing yeast apoptosis network. VCP is important player in can- 
cer cell survival and can be used as a target for cancer therapy. It 
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also serves as detector of aggregated proteins that are known to 
be cause of neurodegenerative diseases such as Parkinson's and 
Alzheimer's. Our model showed that cell survival is mediated 
by VCP via the degradation of IkBa that leads to translocation 
of NFkB to nucleus. This result is in agreement with experi- 
mental findings that VCP is an essential target in the Akt sig- 
naling pathway suggesting that presented model in combination 
with experimental approach would represents a promising plat- 
form to study complex cellular processes involved in cancer and 
neurodegeneration. 

Exploiting the advantages of BN models that enables extrac- 
tion of system level properties of large networks and extensive 
state exploration and converting discrete model to continuous 
model our results suggested that in contrast to human apoptotic 
network yeast apoptosis is linear process whose regulation does 
not involve any complex feedback loops. Analysis of the second 
model showed that with certain adjustments Boolean model of 
yeast apoptosis can be adapted for studies of apoptosis in higher 
organisms. 

Our results show that even without kinetic and qualitative data, 
it is possible to build models that can simulate relevant states of 
yeast physiology and regulations and can contribute to further 
understanding of biology. Since yeast is a preferred model organ- 
ism for many studies of fundamental processes in a Eukaryal cell, 



we argue that in silico studies of yeast will be an important contrib- 
utor to the understanding of complex cellular regulations, such as 
cell death pathways and that the applications will extend toward 
study of regulations or causes of human diseases such as cancer 
and neurodegeneration. 

MATERIALS AND METHODS 
APOPTOSIS NETWORK SETUP 

Based on the extensive literature study and databases search (Sac- 
charomyces Genome Database-SGD; Cherry et al., 201 1 ) the apop- 
tosis network consisting of 73 nodes and 115 edges is constructed. 
These nodes are selected based on their interactions and their sub- 
strates involved in apoptosis. Cell Designer (Funahashi et al., 2003) 
was used to visualize genes and their connection. Cell Designer 
allowed us to represent network using comprehensive graphical 
notation. Moreover, Cell Designer is able to connect to online data 
bases such as KEGG (Kanehisa and Goto, 2000), BioModels (Li 
et al, 2010) and PubMed which expands its connectivity, visual- 
ization and model building. The most prominent advantage of Cell 
Designer is its ability to support SBML (Systems Biology Markup 
Language; Hucka et al., 2003) format. Creating SBML file avoids 
creating logic rule file which is more error prone when rules have 
to be defined by user in simple text and it is not in form of easily 
drawing boxes. 
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FIGURE 8 | Interaction Matrix of Yeast Apoptosis Network. Each row 
corresponds to single species and each column corresponds to reactions. A 
red matrix element e, indicates an inhibition influence of species /on reaction 
j. In contras a green filed shows an activation influence of species / on 
reaction j while a blue box indicates species / gets activated in reaction / and 



as it is expected the black cells indicates that species i does not participate in 
reaction /. Number of interactions where the species is involved is mentioned 
at the end of each row as connectivity number of each species. Number of 
reactions that activates/inhibits is mentioned in brackets (Color coding: 
red - inhibition, green - activation and black - no interaction influence). 
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FIGURE 9 | Dependency Matrix for Yeast Apoptosis Network. In the 

dependency matrix each element m„ represent the relation between an 
affecting and an affected species. Former is specified at the bottom of each 
column and later is shown at the beginning of each row. At intersection of /* 
column and /* row there are three possibilities. A yellow box indicates 
species /is an ambivalent factor meaning both activating and inhibiting path 
exist from species / to species /. Similarly a dark green and or light green cell 
shows a total and a non-total activator respectively. Another possibility is 



COMPUTATION OF LOGICAL STEADY STATES 

Identification of Logical Steady States (LSSs) in a Boolean net- 
work is an important task as they comprise the states in which a 
gene-regulatory network resides most of the time. Strong biolog- 
ical implication can be carried out by LSSs of the network. LSSs 
can even be linked to phenotype (Kauffman, 1969). CNA (Klamt 
et al., 2006) is used to calculate all possible LLSs based on specified 
initial value for each gene and signal flow in network. LSSs are used 
to evaluate the network behavior under perturbation and changes 
in network structure. 

IDENTIFYING NETWORK WIDE DEPENDENCY 

Considering a pair of nodes (a, b) dependency is defined as the 
influence of node a on node b and vice versa, node a influences 
node b as follow: a is a total activator of b if there is an activat- 
ing path between a and b, a is a total inhibitor of b is there is 
an inhibition path from a to b or a is an ambivalent effecter of 
b meaning that there is an intermediate node which is involved 
in a negative feedback loop. Using CNA the interaction matrix 
and dependency matrix were drawn from the network (Figures 8 
and 9). 



having dark or light red indicates species / is a total inhibitor or non-total 
inhibitor of species /. whenever there is no path from spices / to species / 
the intersection cell is field by black. The hyper graph underlying the network 
is a directed graph and consequently the dependency matrix is 
non-symmetric. (Color coding: light and dark green - complete and 
incomplete activation, dark and light red - complete and incomplete 
inhibitor, yellow indicates an ambivalent factor and black indicates that there 
is no dependency between two species). 



CONVERSION OF DISCRETE TO CONTINUOUS MODEL USING SQUAD 

In order to convert the schematic network into Boolean model 
we used SQUAD (Mendoza and Xenarios, 2006), a user friendly 
graphical software which is suitable for modeling signaling net- 
work where kinetic reactions are not available. Simulation in 
SQUAD consists of three steps: (1) the network is first loaded 
from the SBML file. Components of the network are presented as 
nodes and value of each node represents the state of that node. 
SQUAD converts the network to discrete dynamic model. Using 
Boolean algorithms all steady sates in network are calculated. (2) 
network is converted to a continues dynamic model generating 
sets of ordinary differential equations (ODEs) and steady states 
achieved in pervious step. (3) SQUAD allows perturbation in order 
to understand the role of each node within the network. 

Using Reduced Order Binary Decision Diagram (ROBDD), it 
is possible to calculate steady states (Di Cara et al., 2007). ROBDD 
is a memory efficient data structure which is widely used in elec- 
tronic field and has been proven to work for large binary networks. 
Moreover, ROBDD computes steady states for large networks 
(n > 50) in matter of seconds. Another advantage of this algorithm 
is its ability to identify the cyclic steady states. These oscillating 
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states are reachable when system identifies a cyclic pattern instead 
of one single state. 

In the first step, schematic network is supplied to SQUAD 
resulting in a discrete network, generating a set of either cyclic 
or single steady states. These steady states are then used to convert 
discrete model to continues. Given all calculated steady states or 
cycles we examined all possible outcomes for each input. Depend- 
ing on the desire form of output, discrete to continuous conversion 
can be carried out either as a complete or progressive mode. Our 
simulations are performed in the complete mode since cell under- 
going apoptosis should die after certain time and it is expected to 
maintain a constant level of apoptosis or survival (in case when 
apoptosis is not activated) at the end of each run. As an opposite 
to complete mode, the progressive model allows the user to stop 
the simulation at any time even before reaching the steady states. 

CONVERSION OF DISCRETE TO CONTINUOUS MODEL USING ODEFY 

Reconstructing Boolean model of yeast apoptosis from qualitative 
knowledge never gives details about concentration of molecules 



in different time points. For this purpose the discrete Boolean 
model is transformed to continuous model using Odefy. Odefy 
uses the multivariate polynomial interpolation in order to trans- 
form the logical rules into sets of ODEs. Yeast apoptosis Boolean 
model is converted to continuous model using Hill Cube and nor- 
malized Hill Cube where the Hill function is normalized to the 
unit interval. Behavior of biochemical reactions can be seen as 
a sigmoid Hill function represented as /(3c) = x"/(x" + k"). 
Where, n is a Hill coefficient and determines the slope of the 
curve and is a measurement of cooperativity of the interactions, 
and parameter k corresponds to values 1 and 0 in the Boolean 
model in the following manner: threshold value above given k in 
Boolean model is considered as 1 (on state) and below as 0 (off 
state). 
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